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Abstract 

This article deals with flow of plane curves driven by the curvature and external 
force. We make use of such a geometric flow for the purpose of image segmenta- 
tion. A parametric model for evolving curves with uniform and curvature adjusted 
redistribution of grid points will be described and compared. 
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^ ■ 1 Introduction 

> : 

^ . Finding shapes in images and their transformation to some usable form are very impor- 
■ tant and attractive tasks in numerical computation. For such an image segmentation, 
various numerical methods have been proposed. In this article, we will discuss 2D image 
segmentation by using moving plane curves which evolve according to the law: 

normal velocity = curvature -|- external force. (1) 

Under suitable setting of an external force and appropriate choice of numerical parameters, 
segmentation of an edge of the given object image will be obtained as a stable stationary 
solution. As it will be obvious from numerical examples in the following sections, one can 
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obtain a good approximation of the image edge after sufficiently long evolution of planer 
curves. The CPU time is very short, varying from 1 to 10 seconds order. 

Evolving curves can be treated numerically as well as mathematically in several ways: 
for instance, direct approach based methods (Lagrangian method), level-set methods, 
phase-field methods, etc. Our method is based on a direct approach. In the direct ap- 
proach, however, evolving curves are approximated by polygons, the vertices are tracked 
with discrete time. This type of tracking methods often may produce unstable computa- 
tional results. For instance, vertices may accumulate somewhere and may be very sparse 
elsewhere. Hence, a more stable numerical tracking scheme for evolving curves is re- 
quired. To avoid these undesired phenomena, some kind of points redistribution method 
have been extensively studied by many authors. One of the useful ways is redistribution 
of points uniformly along the curve by means of the tangential velocity. In the series of 
papers [H [5l El [7] , K. Mikula and D. Sevcovic proposed asymptotically uniform redistri- 
bution method for a general geometric equation ([T]), and, in particular, they apphed it to 
image segmentation. We also refer to [10] for another approach to image segmentation of 
moving objects, where asymptotically uniform redistribution method has been applied for 
a geometric equation ([1]) with the external force including the probability that the curve 
belongs to the background or target of the segmented sequence of images. For image 
segmentation of indirect approach, we refer, e.g., [HE] for the readers. 

Besides the uniform redistribution method, a curvature adjusted method was proposed 
by D. Sevcovic and S. Yazaki in [11], in which one can find a short history of redistribution 
methods. Their method provides adequate redistribution of grid points depending on the 
absolute value of curvature. The idea comes from the tangential redistribution in the 
case of a crystalline curvature fiow [13j. The advantage of this method is that few points 
are needed to get sharp corners on edges in comparison with the uniform redistribution 
method. 

Organization of the present paper is as follows. In the next section, the evolution 
equations will be introduced. In Section 3, a numerical scheme is proposed and analyzed. 
In the last section, computational experiments are presented. 

2 Evolution equations of curves 

Our method is based on the direct parametric (Lagrangian) approach. According to [llj, 
we introduce the evolution equation as follows. We consider an embedded and closed 
plane curve F which is parametrized counterclockwise by a smooth periodic function 
X : R/Z D [0, 1] R2 such that F = Image(cc) = {x{u); u G [0, 1]} and \dux\ > 0. Here 
and after, the circle M/Z is represented by [0, 1] with periodic condition x{0) = x{l), and 
we define d^F = dF/dC,, and \a\ = ^a.a where a.h is the Euclidean inner product between 
vectors a and 6. The unit tangent vector can be defined as T = duxl\du'x\ = dgX, where 
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s is the arc-length parameter and ds = \dux\du, and the unit inward normal vector is 
defined in such a way that T A N = 1, where a A b = det(a, b) for two dimensional 
column vectors a, b. The signed curvature in the direction N is denoted by k, e.g., the 
sign of k is plus for strictly convex curves. Let u be the angle of T, i.e., T = (cos u, sin u) 
and N = (— sin z/, cos z/). The geometric evolution law ([T]) of plane curves is stated as 
follows: For a given initial curve r° = Image(a;°) = F, find a family of curves {F*}t>o, 
F* = {x{u,t); u G [0,1]} starting from a;(u,0) = x^{u) for u G [0,1] and evolving 
according to the geometric equation 



where v is the normal velocity in the inward normal direction, and F is a given external 
force defined everywhere on the plane. The evolution equation of a position vector x is 
given as: 



Here a is the tangential component of the velocity vector. Notice that the presence of 
an arbitrary tangential component does not affect the shape of curve (see [2l Proposition 
2.4]). Also, dt^ means dtf{u,t), and dsF is given in terms of u: dgf = \duX\~^duF. 

The external force F{x) plays a very important role in the application to the image 
segmentation, because it describes the image to be segmented. In our approach, its value 
is given by the intensity of the image on each pixel. The formula how to calculate the 
value of F from the picture will be mentioned in (fT2l) . 

According to [Ij, equation ([3]) can be written as the following PDE: 



A solution X is subject to the initial condition x{-, 0) = x^{-). The curve is closed, so we 
have to impose periodic boundary conditions of x for u G [0, 1] and z/(l, t) = z/(0, t) + 2tt. 
The term a is responsible for the tangential redistribution. According to [11], [12], the 
curvature adjusted tangential velocity has been proposed in the form: 



V = k + F{x), 



kx = {k + F)N + aT, x{-,0) = x\-). 



(3) 



dtX = dgX + FN + adgX . 



(4) 



(5) 



where L* 



is the curve length in time t, u; is a given positive constant, and 



(p{k) = 1 — £ + £a/1 — £ + e/c^. 



(6) 



/ = ipik)kik + F)- ip\k) [dlk + dlF + k^k + F)) , ^'{k) 



d 

dk 



ip{k), 
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Equation is designed for the exponential convergence of an extended relative local 
length limt^+oo\duX\(f{k)/{L^{(f{k))) = 1, with the exponent —ut. The function f{k) 
given by ([6]) is very important because it controls the redistribution of grid points. Note 
that ip{k) — > 1 if e 0^, and (p{k) — > \k\ if e —>■ The function (p{k) = 1 produces 
the uniform redistribution of points for uj = and asymptotically uniform redistribution 
of points for > (see [H El El Ej), and the function ip{k) = \k\ is used implicitly for 
the crystalline curvature flow p!3]. By choosing e G (0, 1), we obtain curvature adjusted 
redistribution 

To get unique solution a, we assume the following additional condition for a: 

(a(-,t)) = 0. (7) 
The details can be found in the forthcoming paper [T2] . 

3 Numerical scheme 

We follow the numerical scheme developed in [HI [12]. For a given initial A^-sided polygonal 
curve P° = [jf^iS^, let us find a family of A^-sided polygonal curves {'P-'}j=i,2,..., = 
UiLi ^i-: where iS/ = [£C^_i, is the iih. edge with x^^ = x^^ for j = 0, 1, 2, . . .. The initial 
polygon is an approximation of F". We construct as an approximation of F* at 
the jth time t = tj, where to = and tj = J2i=o (j = 1; 2, . . .) with the adaptive time 
increments > for / = 0,1,.... The updated curve {V^^^} is determined from the 
data {V^ at the previous time step by using discretization in space and time of a closed 
system of PDEs (ll) and ([5]). 

The way of discretization is similar to [6], it means that our scheme is semi-implicit 
scheme and the discretization is based on the flowing finite volume approach from [3] . See 
also [8, 9]. Discrete quantities aj, kj, uj, xl, rj = \S-\ ofV^ {i = 1,2, . . . , N) are splitted 
into two categories: kj, and rj take constant value on the ith edge Sj (the finite volume), 
whereas aj and x^ are defined at the ith vertex and take constant values on the zth dual 
edge 5*"' = [a;*"', a;*:J_;^] (the corresponding dual volume) where x*-^ = {x^ + xl_^)/2. The 
duality is defined in such a way that F*-' = (F^ + F^_,_^)/2 (where F = fc, i/, r) take constant 
value on 5*-', and a*-' = {aj + J/2 takes constant value on Sj. Hereafter all quantities 
except {z//} satisfy the periodic boundary condition Fq = Fat, Fat+i = Fi. We will use the 
following notation: ((9^F)^ = (F^"^^ - F^)/r, F^in = mini<i<Ar Fj, |F|max = maxi<i<iv |Fi|, 
(F^) = ^\rl/U, U = ZLri is the total length of P^ {d,f)l = (Fj - fU/rj, 

(dsni = - ^ti)/rl, (ds^m = (fU - H)/rr, ids.f*)l = - H^/rr; and (VO 
means {i = 1,2, . . . , N), whereas (V'i) means {i = 2,3, . . . , N). 

At first, the initial values for } are computed from the initial data ac- 
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cording to 

(i) rf = \Pi\, Pi = {Pn,Pi,) =xl- x\_^ (Vi), 

(ii) = ^ sgn(pi_i A pi+i) arccos j ^'"^^^ ) (^^) 



(iii) ^^ = |™(Pn/rD (P.>0) ^^^^^ 

|^27r - arccos(pj,/r^) (p^a < 0) 



(iv) z/^ = z/i, 



(v) i./ = <''^^^^'' (l^.-^til > l^'.±2vr-z/,til) 
' z/j (otherwise) 



(vi) i^^+i 



z/i±27r (|z/i - z/;^| > |z/i ± 27r - z/^l), 
i/i (otherwise), 



^viij = v{ - ^i/^^i - i/^j. 

Next we discretize equation (JSj) for redistribution term a and compute its values. By 
taking discrete time stepping, we obtain the following closed system of finite difference 
equations of tangential velocity {0^^^}: 

^i^)o^' = ^i^i)o^ + i^l (V'O, (9) 



where 



+ 



r^j y 



where > is a given constant (in the following experiment, uj = 50000 will be used) 
and F*-' = F{x*-'). From we obtain 

i 

^{k:nc4'-' = ^{k;^)at' + ^l ^i = J2^i {Vt). (10) 

1=2 

On the other hand, the equality YliLi "^l^ c>Pi^^ = follows from {a*^^^) = which is a 
discretization of ([7j). Therefore {a{^^}fLi can be determined uniquely from (fTOj) and 



Finally, we discretize equation ([3]) for the position vector x. According to [HI [12], we 
obtain the following closed semi-implicit system: 

{drX)l = {ds.{dsX))l^' + c4'-\ds.X*)l^' + FixDNiu^). 
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This system is equivalent to the following linear system for position vectors xi~^^ subject 
to periodic boundary conditions: 

-ai^Kxtl + (1 + t^T)xl^' - c^Kxi^l =xi + F{xi)N{ur)r (Vz), (11) 
where bl~^^ = a^'^ + c^"*"^, and 

I \ t / I \ t 

The above linear system is a diagonally dominant and solvable if we choose the following 
time step: 

1 4/1 |a^+Vax\ 

V = ^{^ + - ^ , A>0. 



ri(l + A) ' r"' ■ V r-' 

' ^ ' min ^ min 



The transformation formula for mapping pixel colors to F will be given in the next 
section. 

Our scheme is described in the following algorithm. 



1. Set j := 0. Start up with a sufficiently large initial closed curve {P^} given by its 
vertices 

2. Compute {r^}^^, {/c-}^^, {yl^f^^ by means of equations (JED- 

3. Get the external force {F/}^^ and {F^*"*}^]^ from the intensity function lix). 

4. Compute {a^"^^}^^ by solving equation ([2]). 

5. Compute {a;!"^^}^]^ by solving equation ffTTl) . 

6. Set j := j + 1 and go to step 2. 

In our numerical experiment, we always obtain almost stationary shape V\ if j is large 
enough. An appropriate choice of stopping condition for the above algorithm is studied 
in our on-going research. 

4 Experimental comparison between the uniform and curvature adjusted re- 
distribution methods 

In this section, we will describe how to use our method for the image segmentation. 
Moreover, we will compare the uniform and curvature adjusted redistribution of points. 
For our computations, we did not use real images but artificial ones as they are better for 
explaining the properties of curve evolution and points redistribution. We have selected 
three images for this article. The first one (Figure [2]) is a kanji (Chinese character) from 
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Professor Mimura's name. You can see his name in kanji in Figure [H The shape of the 
character 'yasu' is a relatively complex because it contains many gaps and sharp corners. 
The second example image (Figure [7]) contains a lot of noise and the character 'S' is 
damaged by it. Finally, the third image (Figure [9l) contains a silhouette of two people for 
which it is difficult to segment correctly the shape of their faces. 




mi mura masa yasu 



Figure 1: The name of Professor Masayasu Mimura written in kanji. 





-1.5 -1 -0.5 0.5 1 1.5 



Figure 2: The original bitmap image (left) and the very accurate image segmentation at 
t = 0.025. We chose = 1500, e = 0, and F G (-100, 100) (right). 



In all following examples, the initial condition is a circle with the radius 1.5. Then, 
according to the algorithm described in the previous section, we let the curve evolve until 
the stopping condition is reached. It is important to pay attention to proper choice of 
segmentation parameters. 

The image segmentation by this method has three important parameters. The pa- 
rameter N denotes the number of points which the curve consists of. It is clear that the 
more points we use the better final shape we get. On the other hand, higher number 
of grid points makes the computation much slower. The parameter e (see equation (Q) 
controls the curvature adjusted redistribution. For e = we obtain the asymptotically 
uniform redistribution. For e G (0, 1) we obtain a curvature adjusted redistribution. For 
numerical computations, we recommend to choose the value between 0.1 and 0.2. The last 
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but very important function is the external force term F depending on given 600 x 600 
pixels bitmap picture and influences the final shape significantly. 

In the following computations, the target figure is given by a digital gray scale bitmap 
image represented by integer values between and 255 on 600 x 600 pixels. The values 
and 255 correspond to black and white colors, respectively, whereas the values between 
and 255 correspond to gray colors. In our three examples, the target shapes are given in 
white color with black background. 

Given a figure, we can construct its image intensity function I{x) G {0, . . . , 255} C Z 
defined in the computational domain Q := (—1.5, 1.5) x (—1.5, 1.5). We remark that 
I{x) is piecewise constant in each pixel. 

We define the forcing term F{x) as follows: 

-^(•^) Fmax {Fmax Fmin) ^ l-^^) 

where F^ax > corresponds to purely black color (background) and Fmm < corresponds 
to purely white color (the object to be segmented). Maximal and minimal values deter- 
mine the final shape because in general 1/F is equivalent to the minimal radius the curve 
can attain. The choice of small values of F causes the final shape to be rounded or the 
curve can not pass through narrow gaps. Various choices of Fmax and Fmin are shown in 
Figure [5l Figure 5(a) illustrates a situation when the force is too small and the curve can 
not evolve through the obstacle. Increasing the value of F (Figures 5(b) , 5(c) , 5(d) ), the 
shape becomes more and more sharp. 

Now we will describe the results of a uniform {e = 0) and curvature adjusted {e > 0) 
redistribution method. Figures [3] and H] show the evolution of segmentation curves in time. 
In both cases, the number of points = 250, the external forcing term F is calculated 
by (|T2ll where Fmin = —100 and Fmax = 100 but the e is different. One can see that 
for e = (Figure [3]) very sharp corners are smoothed or contain too few points while for 
e = 0.2 (Figure Hj) even the spike under the first horizontal line is kept. Relatively high 
value of the forcing term makes the evolution fast. Both computations were performed 
for the same time interval t G (0,0.025). Figure [2]^ right) shows the result with very high 
number points N = 1500 and e = 0. 

Figure [6] illustrates a comparison of the last step of the evolution for the uniform and 



curvature adjusted redistribution for two different choices of A^. The curve in Figure 6 (a' 



contains 250 grid points. One can clearly see that for e = the narrow gap under the 
second horizontal line is not segmented well and also high curvature parts are rounded. 
For = 350 (Figure 6(b)), the segmentation is quite similar for both cases but still for 
e = 0, sharp parts are not as good as for the case e = 0.2. 

We tested the redistribution method also for images containing noise or some addi- 
tional artifacts. Figure [7| illustrates this case. According to (fT2l) . gray colors do not 
generate so strong force. Hence the curve can easily pass through. If the noise generates 
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(a) t = 0.005 (b) t = 0.010 (c) t = 0.015 (d) t = 0.020 (e) t = 0.025 



Figure 3: = 250, e = 0, F e [-100, 100]. 








(a) t = 0.005 



(b) t = 0.010 



(c) t = 0.015 



(d) t = 0.020 



(e) t = 0.025 



Figure 4: = 250, e = 0.2, F E [-100, 100]. 



too large external force (almost white color) then it prevents the curve to move further. 
But on the other hand, very high curvature will appear and causes the curve to over- 
come the noise. When we work with noisy images, it is very important to choose suitable 
interval for F. If too wide interval is chosen then even a small noise can stop the com- 
putation. On the other hand, too narrow interval will cause bad shape segmentation. In 
our computations for noisy images, we chose F G (—30, 35). 

Since the shape in Figure [7] does not have many details and the length is not so big, 
we do not need to use high number of points A^. The computation in Figure [7t^ right) was 
done for = 800 only for a reference value. The comparison has been performed for 



A^ = 80 (Figure 8(a) and A^ = 150 (Figure 8(b)). 



The last example depicts silhouette of two persons (Figure [9]). There are parts having 
high curvature and the faces of persons have small details. We need a large forcing term 
to find the shape correctly. Since there is no noise, we can do it. Even the simulation for a 



quite low A^ = 150, Figure 10(a) shows that we still get a reasonably good approximation 
with e = 0.15. Computation with e = did not find details near hands and also hat is 
not as sharp as it should be. By increasing A^ to 350, we were able to achieve a good 
approximation with both values of e. But again, sharp edges in the hat area are better 



segmented with e = 0.2 (Figure 10(b)). 
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(a) Fe (-20,30) 



(b) F e (-30,30) 



(c) F e (-40,40) 



(d) F e (-100,100) 



Figure 5: The comparison of several choices of the external force with N = 350 and 
e = 0.2. 

All computations were performed on a standard personal computer with Intel processor 
at 2.4 GHz. The time of computation was never higher than about 10 seconds. The 
computation was faster for £ = than for £ > 0. For N about 200, the CPU time is 
usually less than 1 second. 
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(a) N = 250 (b) N = 350 
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Figure 7: The original bitmap image (left) and a very accurate image segmentation at 
t = 0.1. We chose N = 800, e = 0, and F e (-30, 35) (right). 
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(a) TV = 80 



(b) N = 150 



Figure 8: A comparison between image segmentations with e = and e = 0.2. Here 
F e (-30,35) and t = 0.1. 



12 



Figure 9: The original bitmap image (left) and very accurate image segmentation at 
t = 0.025. We chose = 1500, e = 0, and F G (-100, 100) (right). The original image 
is taken from: http : / /www . zekam . uni-bremen . de/ siluette . gif . 
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(a) TV = 150 (b) N = 350 

Figure 10: A comparison between image segmentations with e = and e = 0.15. Here 
F e (-100,100) and t = 0.05. 
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